# -----------------------------------------------------------------------------#
# title:   Elite-Public Interaction on Twitter:  
#          EU issue Expansion in the Campaign
# content: code to reproduce results from main text
# journal: European Journal of Political Research
#    year: 2020
# -----------------------------------------------------------------------------#
# packages
library("ggplot2")
library("firasans")
library("scales")
library("grid")
library("dplyr")
library("stringr")
library("xtable")
library("brms")

Sys.setenv(TZ = "Europe/Berlin")


set.seed(0101100)

# --------------------------------------------------------------------------#
# a. data loading and added variables
# --------------------------------------------------------------------------#
# load data
load("ejpr-tweets.Rdata")

# additional recoding, helper variables for later modeling
cand_own <- cand_own %>%
  mutate(is_eu_f     = factor(is_eu),
         is_eng_f    = factor(is_eng),
         is_response = ifelse(n_resp_indata == 0, 0, 1),
         politician  = pol_uid)

# --------------------------------------------------------------------------#
# b. multivariate models
# --------------------------------------------------------------------------#
# PLEASE NOTE: the models take quite some time to run

options(mc.cores = 4) # edit as you see fit
# models reported in Table 5
# models without interaction
resp_noint_uk <- brm(n_resp_indata ~ 1 + is_eu_f + is_eng_f + 
                       sitting + electability +
                       fc_2sd + (1 + is_eu_f + is_eng_f | politician),
                     family = "negbinomial", data = filter(cand_own, country == "UK"),
                     chains = 4, cores = 4, iter = 3000)

resp_noint_de <- brm(n_resp_indata ~ 1 + is_eu_f + is_eng_f + 
                       sitting + electability +
                       fc_2sd + (1 + is_eu_f + is_eng_f | politician),
                     family = "negbinomial", data = filter(cand_own, country == "Germany"),
                     chains = 4, cores = 4)

resp_noint_es <- brm(n_resp_indata ~ 1 + is_eu_f + is_eng_f + 
                       sitting + electability +
                       fc_2sd + (1 + is_eu_f + is_eng_f | politician),
                     family = "negbinomial", data = filter(cand_own, country == "Spain"),
                     chains = 4, cores = 4)

resp_noint_gr <- brm(n_resp_indata ~ 1 + is_eu_f + is_eng_f + 
                       sitting + electability +
                       fc_2sd + (1 + is_eu_f + is_eng_f | politician),
                     family = "negbinomial", data = filter(cand_own, country == "Greece"),
                     chains = 4, cores = 4)

# models with interaction
resp_uk <- brm(n_resp_indata ~ 1 + is_eu_f + is_eng_f + is_eng_f*is_eu_f +
                 sitting + electability +
                 fc_2sd + (1 + is_eu_f + is_eng_f | politician),
               family = "negbinomial", data = filter(cand_own, country == "UK"),
               chains = 4, cores = 4, iter = 3000)

resp_de <- brm(n_resp_indata ~ 1 + is_eu_f + is_eng_f + is_eng_f*is_eu_f +
                 sitting + electability +
                 fc_2sd + (1 + is_eu_f + is_eng_f | politician),
               family = "negbinomial", data = filter(cand_own, country == "Germany"),
               chains = 4, cores = 4)

resp_es <- brm(n_resp_indata ~ 1 + is_eu_f + is_eng_f + is_eng_f*is_eu_f +
                 sitting + electability +
                 fc_2sd + (1 + is_eu_f + is_eng_f | politician),
               family = "negbinomial", data = filter(cand_own, country == "Spain"),
               chains = 4, cores = 4, iter = 3000)

resp_gr <- brm(n_resp_indata ~ 1 + is_eu_f + is_eng_f + is_eng_f*is_eu_f +
                 sitting + electability +
                 fc_2sd + (1 + is_eu_f + is_eng_f | politician),
               family = "negbinomial", data = filter(cand_own, country == "Greece"),
               chains = 4, cores = 4)

# same for full sample with country controls

resp_noint_all <- brm(n_resp_indata ~ 1 + is_eu_f + is_eng_f + 
                        sitting + electability + country +
                        fc_2sd + (1 + is_eu_f + is_eng_f | politician),
                      family = "negbinomial", data = cand_own,
                      chains = 4, cores = 4, iter = 3000)

resp_all <- brm(n_resp_indata ~ 1 + is_eu_f + is_eng_f + is_eng_f*is_eu_f +
                  sitting + electability + country +
                  fc_2sd + (1 + is_eu_f + is_eng_f | politician),
                family = "negbinomial", data = cand_own,
                chains = 4, cores = 4, iter = 3000)

# summary(resp_noint_uk)
# summary(resp_uk)
# summary(resp_noint_de)
# summary(resp_de)
# summary(resp_noint_es)
# summary(resp_es)
# summary(resp_noint_gr)
# summary(resp_gr)

# summary(resp_noint_all)
# summary(resp_all)

# model summary in Figure 3 
# themes used in the final version, noted where it can be commented out
all_mat <- as.matrix(resp_all$fit)
# colnames(all_mat)[1:11]
all_diff <- data.frame(rbind(quantile(exp(all_mat[, 1] + all_mat[, 2] +
                                            all_mat[, 3] + all_mat[, 11]) - 
                                        exp(all_mat[, 1] + all_mat[, 2]), c(0.5, 0.025, 0.975)),
                             quantile(exp(all_mat[, 1] +
                                            all_mat[, 3]) - 
                                        exp(all_mat[, 1]), c(0.5, 0.025, 0.975)),
                             quantile(exp(all_mat[, 1] + all_mat[, 2]) - 
                                        exp(all_mat[, 1]), c(0.5, 0.025, 0.975)),
                             quantile(exp(all_mat[, 1] + all_mat[, 2] +
                                            all_mat[, 3] + all_mat[, 11]) - 
                                        exp(all_mat[, 1] + all_mat[, 3]), c(0.5, 0.025, 0.975))))
names(all_diff) <- c("est", "lower", "upper")
all_diff$country <- "Pooled"
all_diff$what    <- c("1. Engaging - Broadcasting",
                      "1. Engaging - Broadcasting",
                      "2. EU - Non-EU",
                      "2. EU - Non-EU")
all_diff$group_on <- c("EU", "Non-EU", "Broadcasting", "Engaging") 



uk_mat  <- as.matrix(resp_uk$fit)
uk_diff <- data.frame(rbind(quantile(exp(uk_mat[, 1] + uk_mat[, 2] +
                                           uk_mat[, 3] + uk_mat[, 8]) - 
                                       exp(uk_mat[, 1] + uk_mat[, 2]), c(0.5, 0.025, 0.975)),
                            quantile(exp(uk_mat[, 1] +
                                           uk_mat[, 3]) - 
                                       exp(uk_mat[, 1]), c(0.5, 0.025, 0.975)),
                            quantile(exp(uk_mat[, 1] + uk_mat[, 2]) - 
                                       exp(uk_mat[, 1]), c(0.5, 0.025, 0.975)),
                            quantile(exp(uk_mat[, 1] + uk_mat[, 2] +
                                           uk_mat[, 3] + uk_mat[, 8]) - 
                                       exp(uk_mat[, 1] + uk_mat[, 3]), c(0.5, 0.025, 0.975))))
names(uk_diff) <- c("est", "lower", "upper")
uk_diff$country <- "UK"
uk_diff$what    <- c("1. Engaging - Broadcasting",
                     "1. Engaging - Broadcasting",
                     "2. EU - Non-EU",
                     "2. EU - Non-EU")
uk_diff$group_on <- c("EU", "Non-EU", "Broadcasting", "Engaging")                   


de_mat  <- as.matrix(resp_de$fit)
de_diff <- data.frame(rbind(quantile(exp(de_mat[, 1] + de_mat[, 2] +
                                           de_mat[, 3] + de_mat[, 8]) - 
                                       exp(de_mat[, 1] + de_mat[, 2]), c(0.5, 0.025, 0.975)),
                            quantile(exp(de_mat[, 1] +
                                           de_mat[, 3]) - 
                                       exp(de_mat[, 1]), c(0.5, 0.025, 0.975)),
                            quantile(exp(de_mat[, 1] + de_mat[, 2]) - 
                                       exp(de_mat[, 1]), c(0.5, 0.025, 0.975)),
                            quantile(exp(de_mat[, 1] + de_mat[, 2] +
                                           de_mat[, 3] + de_mat[, 8]) - 
                                       exp(de_mat[, 1] + de_mat[, 3]), c(0.5, 0.025, 0.975))))
names(de_diff) <- c("est", "lower", "upper")
de_diff$country <- "Germany"
de_diff$what    <- c("1. Engaging - Broadcasting",
                     "1. Engaging - Broadcasting",
                     "2. EU - Non-EU",
                     "2. EU - Non-EU")
de_diff$group_on <- c("EU", "Non-EU", "Broadcasting", "Engaging")       

es_mat  <- as.matrix(resp_es$fit)
es_diff <- data.frame(rbind(quantile(exp(es_mat[, 1] + es_mat[, 2] +
                                           es_mat[, 3] + es_mat[, 8]) - 
                                       exp(es_mat[, 1] + es_mat[, 2]), c(0.5, 0.025, 0.975)),
                            quantile(exp(es_mat[, 1] +
                                           es_mat[, 3]) - 
                                       exp(es_mat[, 1]), c(0.5, 0.025, 0.975)),
                            quantile(exp(es_mat[, 1] + es_mat[, 2]) - 
                                       exp(es_mat[, 1]), c(0.5, 0.025, 0.975)),
                            quantile(exp(es_mat[, 1] + es_mat[, 2] +
                                           es_mat[, 3] + es_mat[, 8]) - 
                                       exp(es_mat[, 1] + es_mat[, 3]), c(0.5, 0.025, 0.975))))
names(es_diff) <- c("est", "lower", "upper")
es_diff$country <- "Spain"
es_diff$what    <- c("1. Engaging - Broadcasting",
                     "1. Engaging - Broadcasting",
                     "2. EU - Non-EU",
                     "2. EU - Non-EU")
es_diff$group_on <- c("EU", "Non-EU", "Broadcasting", "Engaging")  

gr_mat  <- as.matrix(resp_gr$fit)
gr_diff <- data.frame(rbind(quantile(exp(gr_mat[, 1] + gr_mat[, 2] +
                                           gr_mat[, 3] + gr_mat[, 8]) - 
                                       exp(gr_mat[, 1] + gr_mat[, 2]), c(0.5, 0.025, 0.975)),
                            quantile(exp(gr_mat[, 1] +
                                           gr_mat[, 3]) - 
                                       exp(gr_mat[, 1]), c(0.5, 0.025, 0.975)),
                            quantile(exp(gr_mat[, 1] + gr_mat[, 2]) - 
                                       exp(gr_mat[, 1]), c(0.5, 0.025, 0.975)),
                            quantile(exp(gr_mat[, 1] + gr_mat[, 2] +
                                           gr_mat[, 3] + gr_mat[, 8]) - 
                                       exp(gr_mat[, 1] + gr_mat[, 3]), c(0.5, 0.025, 0.975))))
names(gr_diff) <- c("est", "lower", "upper")
gr_diff$country <- "Greece"
gr_diff$what    <- c("1. Engaging - Broadcasting",
                     "1. Engaging - Broadcasting",
                     "2. EU - Non-EU",
                     "2. EU - Non-EU")
gr_diff$group_on <- c("EU", "Non-EU", "Broadcasting", "Engaging")  


all_diff <- bind_rows(all_diff, de_diff, uk_diff, es_diff, gr_diff)

all_diff$country <- factor(all_diff$country, levels = c("Pooled", "Germany", "Greece",
                                                        "Spain", "UK"))


ggplot(filter(all_diff, what == "1. Engaging - Broadcasting"), 
              aes(x = group_on, y = est, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0, linetype = 3, alpha = 0.5) +
  geom_linerange() +
  geom_point(size = 4, shape = 21, fill = "grey60") +
  facet_wrap(~country, ncol = 3, scales = "free") +
  scale_x_discrete("") +
  scale_y_continuous("From broadcasting\nto engaging") +
  scale_fill_discrete("") +
  theme_ipsum_fsc(base_size = 12, grid = "") + # comment out if needed, just a theme
  theme(panel.background = element_rect(fill = "grey98", colour = NA)) +
  theme(strip.text = element_text(size = 16),
        axis.text.y = element_text(size = 14),
        axis.text.x = element_text(size = 12),
        panel.grid.major.y = element_line(colour = "grey90", size = 0.2, linetype = 2),
        panel.grid.major.x = element_line(colour = "grey90", size = 0.2, linetype = 2),
        panel.spacing = unit(1.5, "lines"),
        axis.title.x = element_text(size  = 10, hjust = 0, vjust = 0),
        axis.title.y = element_text(size  = 10, hjust = 1, vjust = 1, angle = 360),
        plot.margin = unit(c(1,1,1,1), "cm"),
        strip.text.y = element_text(angle = 360),
        plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
        legend.position = c(0.9, 0.25),
        strip.background = element_rect(fil = "grey70", colour = "NA"),
        strip.text.x = element_text(hjust = 0.5, colour = "white")) +
  labs(caption = "Note: Median difference between predicted counts with 95% credible intervals as line ranges. 
       Please note the varying y-axis scales used for easier readability.")


# models reported in Table 7
eng_uk <- brm(is_eng ~ 1 + is_eu_f + 
                sitting + electability +
                fc_2sd + (1 + is_eu_f | politician),
              family = bernoulli(), data = filter(cand_own, country == "UK"),
              chains = 4, cores = 4)
# summary(eng_uk)
eng_de <- brm(is_eng ~ 1 + is_eu_f + 
                sitting + electability +
                fc_2sd + (1 + is_eu_f | politician),
              family = bernoulli(), data = filter(cand_own, country == "Germany"),
              chains = 4, cores = 4)
# summary(eng_de)
eng_es <- brm(is_eng ~ 1 + is_eu_f + 
                sitting + electability +
                fc_2sd + (1 + is_eu_f | politician),
              family = bernoulli(), data = filter(cand_own, country == "Spain"),
              chains = 4, cores = 4, iter = 3000)
# summary(eng_es)
eng_gr <- brm(is_eng ~ 1 + is_eu_f + 
                sitting + electability +
                fc_2sd + (1 + is_eu_f | politician),
              family = bernoulli(), data = filter(cand_own, country == "Greece"),
              chains = 4, cores = 4)
# summary(eng_gr)
eng_all <- brm(is_eng ~ 1 + is_eu_f + 
                sitting + electability + country + 
                fc_2sd + (1 + is_eu_f | politician),
              family = bernoulli(), data = cand_own,
              chains = 4, cores = 4)
# summary(eng_all)

# export all (uncomment, if desired)
# save(cand_own,
#      eng_uk, eng_de, eng_es, eng_gr, eng_all,
#      resp_uk, resp_de, resp_es, resp_gr, resp_all,
#      resp_noint_uk, resp_noint_de, resp_noint_es, resp_noint_gr, resp_noint_all,
#      file = "models-fitted.RData")


# --------------------------------------------------------------------------#
# c. descriptive tables in text
# --------------------------------------------------------------------------#
# please note: Figure 1 is generated based on a different dataset and script 
#              see ejpr-fig-1.R
#
#              Figure 2 is not data based figure (no code supplied)

# tables 1 and 3
cand_own %>% group_by(country) %>% 
  summarize(n_polit  = length(unique(politician)),
            n_part   = length(unique(party_nest_id)),
            tweets   = n(),
            is_eu    = mean(is_eu),
            is_eng   = mean(is_eng),
            resp     = mean(is_response),
            rest_total = sum(n_resp_indata))

# table 4
cand_own %>% group_by(country, is_eng) %>% 
  summarize(eu = sum(is_eu == 1),
            neu = sum(is_eu == 0),
            resp_eu = sum(n_resp_indata[is_eu == 1]),
            resp_neu = sum(n_resp_indata[is_eu == 0]),
            ratio1 = resp_eu/eu,
            ratio2 = resp_neu/neu)

# table 6
cand_own %>% group_by(country, is_eu) %>% 
  summarize(eng = mean(is_eng))
cand_own %>% group_by(is_eu) %>% 
  summarize(eng = mean(is_eng))

# table 8
cand_own %>% group_by(country) %>% 
  summarize(all_rt  = sum(n_retweet_public),
            mean_rt = mean(n_retweet_public),
            mean_rt_eu    = mean(n_retweet_public[is_eu == 1]),
            mean_rt_pol   = mean(n_retweet_public[is_eu == 0]),
            mean_rt_broad = mean(n_retweet_public[is_eng == 0]),
            mean_rt_eng   = mean(n_retweet_public[is_eng == 1]))

# --------------------------------------------------------------------------#
# Appendix. sessionInfo
# --------------------------------------------------------------------------#
# R version 3.5.2 (2018-12-20)
# Platform: x86_64-apple-darwin15.6.0 (64-bit)
# Running under: macOS  10.15.3
# 
# Matrix products: default
# BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
# LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
# 
# locale:
#   [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
# 
# attached base packages:
#   [1] grid      stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
#   [1] brms_2.10.0      Rcpp_1.0.3       xtable_1.8-4     stringr_1.4.0    dplyr_0.8.3     
# [6] scales_1.1.0     firasans_0.1.0   hrbrthemes_0.7.2 ggplot2_3.2.1   
# 
# loaded via a namespace (and not attached):
#   [1] Brobdingnag_1.2-6    gtools_3.8.1         StanHeaders_2.19.0   threejs_0.3.1       
# [5] shiny_1.4.0          assertthat_0.2.1     stats4_3.5.2         yaml_2.2.0          
# [9] gdtools_0.2.1        backports_1.1.5      Rttf2pt1_1.3.7       pillar_1.4.2        
# [13] lattice_0.20-38      glue_1.3.1           extrafontdb_1.0      digest_0.6.22       
# [17] promises_1.1.0       colorspace_1.4-1     htmltools_0.4.0      httpuv_1.5.2        
# [21] Matrix_1.2-15        plyr_1.8.4           dygraphs_1.1.1.6     pkgconfig_2.0.3     
# [25] rstan_2.19.2         purrr_0.3.3          processx_3.4.1       later_1.0.0         
# [29] tibble_2.1.3         bayesplot_1.7.0      DT_0.10              shinyjs_1.0         
# [33] withr_2.1.2          lazyeval_0.2.2       cli_1.1.0            magrittr_1.5        
# [37] crayon_1.3.4         mime_0.7             ps_1.3.0             evaluate_0.14       
# [41] nlme_3.1-137         xts_0.11-2           pkgbuild_1.0.6       colourpicker_1.0    
# [45] prettyunits_1.0.2    rsconnect_0.8.15     tools_3.5.2          loo_2.1.0           
# [49] lifecycle_0.1.0      matrixStats_0.55.0   munsell_0.5.0        callr_3.3.2         
# [53] compiler_3.5.2       systemfonts_0.1.1    rlang_0.4.1          ggridges_0.5.1      
# [57] rstudioapi_0.10      htmlwidgets_1.5.1    crosstalk_1.0.0      igraph_1.2.4.1      
# [61] miniUI_0.1.1.1       base64enc_0.1-3      rmarkdown_1.16       gtable_0.3.0        
# [65] inline_0.3.15        abind_1.4-5          markdown_1.1         reshape2_1.4.3      
# [69] R6_2.4.1             gridExtra_2.3        rstantools_2.0.0     zoo_1.8-6           
# [73] knitr_1.25           bridgesampling_0.7-2 fastmap_1.0.1        extrafont_0.17      
# [77] shinythemes_1.1.2    shinystan_2.5.0      stringi_1.4.3        parallel_3.5.2      
# [81] tidyselect_0.2.5     xfun_0.11            coda_0.19-3 